Bounding biomass in the Fisher equation 

Daniel A. BirchJ*] Yue-Kin Tsang[] and William R. Young * 
Scripps Institution of Oceanography, 
University of California at San Diego, 
La Jolla, CA 92093-0213, USA 
(Dated: February 4, 2008) 

Abstract 

The FKPP equation with a variable growth rate and advection by an incompressible velocity 
field is considered as a model for plankton dispersed by ocean currents. If the average growth rate 
is negative then the model has a survival-extinction transition; the location of this transition in 
the parameter space is constrained using variational arguments and delimited by simulations. The 
statistical steady state reached when the system is in the survival region of parameter space is 
characterized by integral constraints and upper and lower bounds on the biomass and productivity 
that follow from variational arguments and direct inequalities. In the limit of zero-decorrelation 
time the velocity field is shown to act as Fickian diffusion with an eddy diffusivity much larger 
than the molecular diffusivity and this allows a one-dimensional model to predict the biomass, 
productivity and extinction transitions. All results are illustrated with a simple growth and stirring 
model. 
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I. INTRODUCTION 

Fisher [l| and Kolmogorov, Petrovskii and Piskunov [2] introduced a partial differential 
equation, now called the FKPP equation, modelling the spread via diffusion of an advan- 
tageous gene through a dispersed population. Skellam \^ applied the FKPP equation to 
the change in abundance of organisms in space and time. Oceanographic applications, par- 
ticularly the dynamics of plankton populations, motivate extending the FKPP model by 
inclusion of an incompressible velocity u(x,t). Thus the FKPP equation considered here is 

P t + u-VP = -fP-r]P 2 + kV 2 P, (1) 

where P(x, t) is the concentration of phytoplankton. This is the simplest model containing 
the four essential ingredients of advection, growth, saturation and diffusion. Because of 
environmental variability, the growth rate j(x,t) may depend on both location x and time 
t. The small-scale diffusivity k and the saturation coefficient rj are taken to be positive 
constants. Figure [1] shows a snapshot of a typical solution of (CQ). 

If 7 in (JT]) is a positive constant then a small initial population grows to occupy the 
entire domain so that ultimately the concentration is uniform (i.e. lim^ooP = 7/77). The 
most interesting aspect of this special case is the interaction between frontpropagation and 
advection which occurs on the way to the uniform steady state [e.g. U, |5j. On the other 
hand, in Figured! u(x,t) continually stirs the population relative to the spatially variable 
growth rate and carrying capacity, resulting in a non-trivial statistically steady solution 

In oceanography a great deal of effort has gone into studying "plankton patchiness" 
and the small-scale structures produced by lateral stirring [e.g. [7|. Here we avoid the issue 
of defining a patch or patchiness, and we also largely avoid considerations of the small-scale 
structure evident in Figure [U Instead we consider a more basic question: can we predict or 
constrain the total amount of plankton in solutions of ([1])? To pursue this goal we develop 
upper and lower bounds on the plankton biomass which depend only on the gross properties 
of the growth rate and stirring. In addition to bounds on the biomass, we also develop 
bounds on the productivity, which may be related to the variance of the concentration. 
These bounds are obtained using mathematical techniques with parallel applications to the 
Navier-Stokes equation and the passive scalar problem [e.g. 

Hoi- 

In Section [U we describe the model growth rate and velocity used to illustrate our general 
results and we make a comparison between an inert scalar and the reactive tracer P in (CQ). 
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FIG. 1: Center: A snapshot of ln(P) from a simulation of (pQ). The plankton concentration along 
the line x = is shown on the left and the growth rate 7(y) is on the right. This simulation uses 
the model in section ITT1 with T = 0.1, £7* = 10, m = 1, r* = 0.25 and k, = lx 10~ 4 . 

The survival-extinction transition is discussed in Section IIIII In Sections IIVI and [V] we 
develop inequalities which constrain the biomass and productivity. Section IVII contains the 
conclusions and discussion. Some mathematical details are contained in three appendices. 



II. AN ILLUSTRATIVE GROWTH RATE AND VELOCITY FIELD 

Our main results will apply to a variety of space- and time-dependent growth rates and 
flows. However, for the sake of illustration, we will present examples using a model defined 
on a doubly periodic square domain, with x, y G [— vr£, ir£), and the growth rate given by 
the "sinusoidal" model: 

l{y) = 7 max P + (1 - T) cos (y/£)] , ye l-n£,ne). (2) 

T G (— oo, 1] is non-dimensional and controls both the average and the spatial structure of 
the growth rate as shown in Figure El 
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FIG. 2: The model growth rate ([2]) for various T. 
Stir ring is provided by a popular model of a random two-dimensional velocity field [e.g 



ring is prov 



12L Il3l . 1 14 , Il5l |. The velocity alternates between t> = and u = 0: 



>/2[/ (cos(k m y + 4> x ) , 0) , for nr < i < (n + l/2)r ; 
u(x,t)={ (3) 

\/2C/(0, cos(A; m :r + 0J) , for (n+l/2)r <t < (n+l)r. 

Above, fc m = m/£ where m is an integer controlling the scale separation between the domain 
scale £ and the velocity field. The phases (fi x and 4> y are randomly chosen each period with 
uniform density on (0,27r). 

The average squared velocity components of the stirring model ([3]) are (u 2 ) = (v 2 ) = \U 2 ■, 
and the flow is homogeneous and isotropic in the sense that 

(uiUj) = \U 2 5ij . (4) 



The angled braces in (jl]) indicate a space-time average and are explicitly defined in (fl9l) . 
Since the flow is isotropic, the single-particle eddy diffusivity may be found by using the 
relationship 

2Dr = ((Ax) 2 ) = ((Ay) 2 ) , (5) 

where Ax and Ay are the x- and ^-displacements of a particle during the time interval t = 
to r. Thus the eddy diffusivity of (J2J) is 

U 2 T 

• (6) 

In addition to the eddy diffusion of individual particles, we also consider the stretching 
and compression of an infinitesimal material line element The length of the element, 



grows exponentially at a rate estimated by 

A = lim £ -1 E [In |£|] , (7) 



where A is the Lyapunov exponent of the flow in ([3]) and E denotes the expectation obtained 
by averaging over an ensemble of material elements. Dimensional considerations show that 
for ([3]) the Lyapunov exponent A has the form 

A = Uk m A(r u ) , r u = Uk m r . (8) 

The non-dimensional parameter t u is the ratio of the de-correlation time r to the shear 
(£7/c m ) _1 . An estimate of the the non-dimensional function A(r u ), obtained by the Monte- 
Carlo method summarized in appendix [XJ is shown as the solid curve in Figure [3](a). The 
dashed curve in Figure [3](a) is the approximation: 

This particular functional form is suggested by analytic solution of closely related problems 

H- 

In addition to m, T and t u , the model is controlled by two more parameters: the Peclet and 
Damkohler numbers. Scaling length with I and time with 7~a X > the Peclet number emerges 
as the ratio of the diffusive time scale £ 2 /k to the advective time scale E/U: Pe = £U/k. The 
Damkohler number is the ratio of the advective time scale i/U to the reactive (or biological) 
time scale 7 m a X : Da = i^m&x/U. For convenience, we use the equivalent parameters 

k* = Da^Pe" 1 = — — , U* = Da" 1 = . (10) 

The non-dimensional renovation cycle length is denoted r* and may be calculated from the 
other non-dimensional parameters: 

r^ 7max r = ^. (11) 

Finally, the saturation constant rj may be completely removed from all equations by scaling 
P with 7 max /?7. 

Using the split-step lattice method of Pierrehumbert [14j , one can efficiently solve both 
(CO) and the forced inert passive scalar (inert scalar from now on) equation 

C t + u.VC = kV 2 C + cos k x y, (12) 




FIG. 3: (a) The solid curve shows a Monte-Carlo estimate of the non-dimensional Lyapunov 
function A(t u ) defined in ([8]); the dashed curve is the fit in ([9]). (b) Spectra of the inert scalar 
C in (|12|) and of the biological tracer P obtained from (pQ). The spectra are "compensated" by 
multiplying the variance spectrum by the wave number so that a k^ 1 Batchelor spectrum appears 
flat. The plankton are growing according to the sinusoidal growth rate in ([2]) with T = 0. The 
P-spectra are line spectra at y/i = 0, where 7 = 7 max , and at yjt = n/2, where 7 = 0. C is well 
described by the Batchelor spectrum, while P has a steeper spectral slope. Panels (c) and (d) show 
snapshots of the P and C respectively; note in panel (c) the contour interval is logarithmic (i.e., 
panel (c) shows Z = ln(P)). Both simulations use the parameters k* = 10 -7 , [/* = 1, r* = 2.2214 
and m = 1. 

with the velocity field in ([3]). Snapshots of simulations with m = 1 in ([3]) and a resolution of 
4096 x 4096 are shown in the bottom panels of Figure [3j The inert scalar C in panel (d) has 
a classic fc _1 -spectrum as shown in panel (b). With the sinusoidal growth rate in ([2]), 
the statistics of the biological tracer P are spatially inhomogeneous and so in Figure E^b) 
we show one- dimensional P-spectra obtained along the lines at which 7 = 7 max and 7 = 0. 
The P-spectra have slopes of —1.26 and —1.40 for the transects with 7 = 7 max and 7 = 0, 



6 



respectively. The inert scalar 's spectral slope is —1.01. All slopes were calculated over the 
range of wave numbers k — 1 to 32. 

Deviations of the P-spectra from the spectra of inert scalars in the ocean are generally 
interpreted to mean that biological processes such as growth and grazing are strongly affect- 
ing the plankton distribution. Continuous sampling methods have allowed comparisons of 



the spectra of chlorophyll in the ocean to the spectra of physical properties since 1972 18], 



but unfortunately there is still no consensus on how the spectra should vary [6( . However, 



a recent paper [19j presents strong evidence that stirring dominates the distribution of phy- 
toplankton along isopycnals on intermediate scales (10 - 100 km). An added complication 
in the ocean is that the "passive" tracer most often compared to chlorophyll is temperature, 
which is dynamically active. 

III. THE SURVIVAL EXTINCTION TRANSITION AND STRANGE EIGEN- 
FUNCTIONS 

Before developing bounds on the plankton biomass and productivity, we first consider 
a more fundamental question: do the plankton survive at all? Early work on this issue 



mc 

a 



udes classic papers on the problem of "critical patch size" in models without advection 
2o| . More recent work ll|, 21, 22, 23] applies to models with advection. We follow 
ll| by considering flows such as ([3]) with non-zero Lyapunov exponent 



Bayly's approach 
and linearizing (Tj[|): 

P t + u-VP = 7P + kV 2 P . (13) 

If the initial plankton concentration is very low everywhere then the quadratic nonlinearity 
is negligible and (I13p governs the initial behavior of the system. 

With the support of numerical simulations (see Figure H]), we assume that following a 
transient stage the evolution of P takes the form 

P(x,t) = e st P(x,t). (14) 

In ( fl4l) the spatial distribution of P is described by P(x,t), the statistically stationary 
"strange eigenmode." The amplitude of the solution either grows or decays according to the 
sign of the "survival exponent" s [ill. \l3\. 
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FIG. 4: (a) The instantaneous biomass P(t) shows the exponential growth of the strange eigen- 
mode in (I14p followed by the dynamic equilibrium when the non-linear saturation term becomes 
important. The initial condition is P(x,0) = 10~ 7 and the non-dimensional parameters are noted 



at the bottom of the figure. The dashed line is the theoretical maximum P(t) = i-'(0)e 7max * 
(b) Increasing [/* in ([3]) decreases s* and results in s* < (extinction) if [7* is greater than about 
6. 

If P is initially positive everywhere, then it will remain so and we can define 

Z = \nP. (15) 

This is a crucial difference between P and the much-studied problem of the decay of concen- 
tration anomalies of an inert scalar to zero. Factoring P into an amplitude e s * and a strange 
eigenmode P(x,t) as in <HM is equivalent to writing 

Z(x,t) = st + Z(x,t), (16) 

where Z = In P is statistically stationary. In terms of Z the problem (jl~3]) is 

Z t + u ■ VZ = 7 (aj, t) + K (V 2 Z + \VZ\ 2 ) . (17) 
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We begin our analysis of (fI7I) by introducing a spatial average denoted by an overbar 
and denned by 

f(t) = 4- I f(x,t)dx, (18) 
where the integral above is over the domain Q with area Aq. For statistically stationary 
fields, such as Z(x,t), we also employ a space-time average defined by 

(/) = lim —J— ! [ f(x,t)dxdt. 
Substituting ffl6l) into ffTTl) and spatially averaging yields 



(19) 



s + Z t — 7 + k 



VZ 



(20) 



Time averaging (I20l . we obtain a fundamental connection between the survival exponent s, 
the average growth rate (7) and VZ: 



{j) + K 



VZ 



(21) 



Thus if (7) > then s > and the population survives (see also 11]). However the converse 
is not true: because of the final term in ( 12TT) . the survival exponent s can be positive even 
if (7) < 0. In Section Till Al we go beyond [11] and explore survival when (7) < 0. 



An adverse environment 



the role of diffusion 



From a theoretical point of view, (7) < is the interesting case: the population might 
survive even though the average environment is adverse. This is illustrated in Figure H^a) 
where the population survives with T = — 1. The role of the diffusive term k (| VZ| 2 ) in (I2T]) 
is quite confusing in this case and the variation of s with k depends on the details of the 
flow. In Figure H^a) decreasing k increases s. 

The limit of large diffusion and consequent extinction is straightforward: if k is very large 
then the population rapidly diffuses over the entire domain and the negative average growth 
rate prevails so that s < 0. In fact, a simple perturbation expansion around n' 1 = quickly 
shows that 



lim s 



(7) 



and 



lim k 



VZ 



0. 



(22) 



As always, the other limit, k — > 0, is potentially singular and holds the possibility that s > 
because the final term in (I2T!) is non-zero. We investigate this possibility by consideration 
of some special cases and via a variational approach. 




FIG. 5: (a) The first three eigenvalues of (123f) . Extinction occurs when si < 0, which for = 0.1 
means V < —20. (b) The solid curve is the f7* = extinction transition computed by numerically 
solving the linearized stability problem ([23]) with s* = 0; the dashed curve is the approximation in 



B. The case (7) < and u = 

As an elementary illustration of survival in an adverse environment with k — > 0, consider 
(I13p with u = 0. In this case (1131) has non-strange eigensolutions, determined by requiring 
P t = and substituting (fl4|) into (ITB1 . The resulting Sturm- Liouville eigenproblem is a 
form of Mathieu's equation, which we express in non-dimensional variables: 

[K*dy + T + (1 — T) cosy] P — s*P , (23) 

where = s/7 max . Because the differential operator on the left of (1231) is self-adjoint, all 
of the eigenvalues s n are real. Figure 0(a) shows the first three eigenvalues as a function of 
— T for the case /t* = 0.1. 

The extinction transition occurs when the largest eigenvalue s\ passes through zero. 
Therefore, to determine the extinction transition systematically we set = in (1231) and 
regard as a new eigenvalue. The largest eigen-K* is the critical value of k*, above which 
extinction occurs. The solid curve in Figure [5](b) marks the boundary between survival 
and extinction in (— F, K*) parameter space computed in this fashion. The dashed curve in 
Figure E^b) is an approximation to the extinction transition obtained using perturbation 
theory in Appendix [Bj 

2 + 6|F| 

* 4|r| + 3r2 • [Z ^> 
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Now we examine the integral constraint (I2ip in light of this example. If — > with T 
fixed and negative then the system enters the survival region in Figure [5](b) . Thus in this 
limit the term ^|VZ| 2 ^ is both non-zero and crucial in ensuring that > 0. Notice 
that since = 0, this singular — ► limit does not involve gradient amplification by 
exponential stretching. 

Instead, we understand the — > limit by finding an approximation to P and directly 
evaluating (^\'VZ\ 2 ^. The first step is calculating the largest eigenvalue of ff23|) using the 
results in Appendix [B] with K = kJ(1 - T) and E = (s* - T)/(l - T) in flB5|) : 



J-(8l _^Mi_^). (25) 

As k» the growth rate of the mode approaches the maximum of 7, namely s* — > 1 



11]. Figure [6]( a) shows that (|25|) is a very good approximation to s* over a wide range of 
and improves as — > 0. The Gaussian approximation of Appendix (EJ used to obtain 
(J2SJ) , is shown in Figure MJo) and is an excellent approximation to P in the region where the 
eigenmode is concentrated. However, to reconcile — > 1 with the integral constraint (T2"Tj) 
we need an approximation which is valid where P y /P is large; ironically, this is the region 
where 7(7/) — s < and P is very small. We can use the method of Wentzel, Kramers and 
Brillouin (WKB hereafter) to obtain the required approximation by assuming that ~ 1 



and re-casting (|23|) in Schrodinger form 24| : 



K*P yy = R(y)P, i? = (l-r)(l-cosy). (26) 
The WKB solution to (1261) which is symmetric about y = tt is 

Avkb = CR'^(y) cosh (V/ 2 • (27) 

Now we can evaluate the term ^|VZ| 2 ^ in fT2"TT) and find that it is independent of ft*: 

/ / - \ 2\ 

ivzi 2 ) = «» P ^ 




tanh 2 (^^V^w) dy> + 0{k\' 2 ) 



« l-r + O^ 2 ), (28) 

where the tanh 2 can be replaced by 1 to evaluate the integral with errors of order k 1 J 2 . 
Equation (1281) reconciles (1251) with (l2Tj) in the singular limit — > 0. Figure [6](c) shows that 
( 127|) provides good approximation to over the entire domain. 
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FIG. 6: (a) The eigenvalue s* in (|23|) as a function of K*. The approximation (|25p improves as 
— > 0. (b) The eigenmode -P(y) in (|23p and the Gaussian approximation from Appendix [Bj The 
WKB approximation ()27[) has a singularity near the origin and is not shown, (c) Three solutions 
for Py/P- Both the Gaussian and WKB approximations are good near the origin, but the Gaussian 
approximation fails in the region where 7(2/) — s < and is therefore not useful for evaluating the 
average in (f2Tj) . 

C. The survival extinction transition in the limit of rapid decorrelation 

In the limit of a rapidly decorrelating velocity the results of the previous section can 
be adapted to make another quantitative prediction of the survival-extinction transition. 
This rapid-decorrelation limit is achieved by taking r* — ► and U* —>■ 00 so that the non- 
dimensional version of the eddy diffusivity in (JS)), namely = U%t*/8, is fixed. Provided 
that 

£>*>/«*, (29) 
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FIG. 7: (a) The biomass (P) as a function of £/* as determined by Monte Carlo simulation. The 
extinction transition is at around U* = 7.2. (b) An expanded view of the data in panel (a), (c) 
Summary of a suite of simulations in which D* = [7^r*/8 is fixed as r* is varied. The survival- 
extinction transition is located via repeated simulations with varying £/*, as in panels (a) and (b). 
As t* is decreased the survival-extinction transition approaches the solid curve previously shown 
in Figure [5jb) with playing the role of . The velocity field is ([3]) with m = 8 to ensure scale 
separation. 

and that there is scale separation between the velocity field and the domain 

m > 1, (30) 

it is plausible that the earlier eigensolution fl24|) . with replaced by in fl23l) . can be used 
to locate the extinction-survival transition in the (— r, D*) parameter plane. Numerical 
simulations, summarized in Figure [TJ show that this eddy-diffusion closure works well. 
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D. Prohibition of extinction: a lower bound on s 



Aside from the rapid-decorrelation limit we do not have a simple means of determining 
the location of the extinction transition in the parameter space. However, in this section we 
use a simple variational method to locate a region of the (k, (it 2 ))-plane where extinction 
is impossible, even if the average growth rate is negative. The region we find is contained 
within the (possibly larger) actual survival region. In other words, we obtain a sufficient 
condition for survival which applies to any incompressible velocity field without restriction 
to rapid decorrelation or scale separation. 

We begin by substituting (fTE]) into (IT7|) and then multiplying by h 2 (x), where h(x) is an 
arbitrary real function of x (but not t). After space-time averaging we have: 



h 2 ) = (/i 2 7 ) - (h 2 u ■ VZ) - nlVh 2 ■ VZ) + K{h 



vz 



(31) 



The first important consequence of insisting that h(x) is independent of t is that ( h 2 Z t ) = 0. 
Rearranging (IHTj) yields a quadratic in h'VZ, and so we complete the square: 



(h 2 ) = (h 2 1 )+K 



hVZ -Vh- 



hu 

2k 



2 


_, hu 




V/i + — 




2k 



(32) 



Dropping the square term containing VZ from the right-hand side of (132!) results in the 
inequality 

s (h 2 ) > </i 2 7 ) - K (| V/i| 2 ) - Ik- 1 (h 2 \u\ 2 ) . (33) 

Survival (s > 0) is guaranteed if we can find any real function h{x) which makes the right- 
hand side of fl33|) positive. 

To avoid guessing at h(x) we apply variational calculus to (133]) and find the function 
h(x) which maximizes the right-hand side. Thus we apply the constraint 



<" 2 > = i 

with a Lagrange multiplier s and maximize the functional 

?[h] EE </l 2 7 > - K <| Vhf) - \K X (h 2 \U\ 2 ) - S (h 2 - 1) 



(34) 



(35) 



The second important consequence of taking h(x) independent of t is that the time-average 
in (h 2 \u\ 2 ) applies only to \u\ 2 , so that for statistically homogeneous and isotropic flows 
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FIG. 8: (a) A schematic illustration; below the dashed curve the sufficient condition (|38p guarantees 
survival. However, the actual survival region may be much larger, as indicated by the region below 
the solid curve, (b) Below each solid curve the sufficient condition (|38p . obtained by numerical 
solution of (|37p . guarantees survival for the indicated value of T in the sinusoidal growth model 

(SHIF), such as the model in ©, (h 2 \u\ 2 ) = (\u\ 2 ) ih 2 ) = U 2 (h 2 ). Thus in this case 

J[h} = {h 2 1 )-K{\Vh\ 2 )-\U 2 K- 1 {h 2 )-~s{h 2 -l) , (for SHIF). (36) 

The corresponding Euler-Lagrange equation is: 

f=0, kVU = (37) 

where h(x) is the function which optimizes (13"3"j) by maximizing the right hand side. 
Conveniently, if we multiply (1371) by h and average, we obtain 

and therefore s > s , (38) 



1 



h 



i.e., the population survives if the maximum eigenvalue s of (1371) is positive. Note that (J3j 
is the same eigenproblem solved in (|23|) to determine the survival-extinction transition with 
u = 0: the effect of u ^ is the same as reducing the growth rate 7(1/) by U 2 /4k. 

The region below the dashed curve in Figure (Hl(a) is where (|38l) forbids extinction; this 
region is contained within the actual survival region for our growth and stirring model, 
which is the area below the solid curve. Note that the survival region associated with the 
solid curve in Figure [8](a) abuts the U 2 -axis, which is qualitatively different from the region 
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below the dashed curve (see also Figure E^b)). The schematic solid curve in Figure [8](a) is 
based on our experience with simulations such as those in Figure H] indicating that if U* is 
below some threshold, roughly Z7* = 6 in Figure IU(b), then the population survives in the 
limit — ► with U* fixed. This behavior is indicated in Figure H^a): decreasing with 
fixed U* = 2 increases s towards Bayly's upper bound 7 max - 

The discrepancy between the actual survival region and the sufficient condition obtained 
from (1331) depends on the details of the flow and growth rate. For example, using the growth 
rate and flow ([3]), the actual survival region is much larger than that calculated from 
( 1331) . This is shown in Figure IU(a) where the plankton survive with T = -1, Ul = 4 and 
0.01 < k* < 1; these parameter values are far above the curve T = — 1 in Figure MJo). Thus 
unfortunately (13"3"j) is not in general a tight lower bound on the survival exponent s. 

IV. THE STATISTICAL STEADY STATE 

We now suppose that the population survives and turn to the statistical steady state which 
ensues once the quadratic nonlinearity in ([1]) halts the exponential growth of the strange 
eigenmode. Two descriptors of this equilibrium are the biomass 23 and the productivity J 3 
defined by 



where () is the space-time average defined in ( jT9l) . An important integral constraint is 
obtained by averaging flTJ): 



Equation (HOI) has the obvious interpretation that in statistical equilibrium reproduction is 
balanced by mortality. Using (1401) we see that the variance of P(x,t), that is (P 2 ) — (P) 2 , 
is equal to ri~ lr P — 23 2 . Thus 23 and CP provide the mean and variance of P(x,t). This is a 
strong motivation for regarding 23 and T as the most fundamental statistical descriptors of 
the system and for attempting to understand their dependence on k and the properties of 
u{x, t) and j(x, t). 

We obtain a second integral constraint by making the change of variables Z = In P in 



03 = (P) , and 7 ee ( 7 P) , 



(39) 




(40) 



Z t + u- VZ = j(x,t) -7]P + k(V 2 Z +\VZ\ 2 ) . 



(41) 
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Averaging (14~TI) gives the equilibrium analog of (I2TT) : 

t 1 (P) = { 1 ) + k{\VZ\ 2 ) . (42) 

This shows that the H> is always greater than 7]~ l (7). This lower bound is only useful if 

(7) > 0, so 

rf£> > max{0, (7)} . (43) 

To obtain a lower bound on the productivity we combine the Cauchy-Schwarz inequality, 
(P 2 ) > (-P) 2 ; with the definition of "? and the identity ( 1401) to obtain (7) x max{0, (7)} < 
We can also employ Cauchy-Schwarz to find an upper bound on the productivity: 
y = (jP) < a/ (7 2 ) (P 2 )- Using fT40|) to replace (P 2 ) by r/^T, and squaring the resulting 
inequality, we obtain the upper bound rf$ < {'J 2 ). Thus, to summarize: 

(7) x max {0, (7)} < ip < <7 2 ) . (44) 

The simple bounds in f l43|) and ( 1441) involve neither u(x,t) nor k. In the next section we 
obtain a more elaborate bound which depends on u(x, t) and k and which applies to the 

case (7) < 0. 



A. A second lower bound on the biomass 

To obtain a lower bound on B e (P) we follow the calculation in subsection UIIIA 
Multiplying (|4T|) by h 2 (x) and space-time averaging we obtain a steady state analog of (|3~T1) : 

v (h 2 P) = </i 2 7 ) - (h 2 u ■ VZ) - k {Vh 2 ■ VZ) + k (h 2 \ VZ\ 2 ) . (45) 

Repeating the manipulations in fl32|) and fl33l) . and again assuming that u(x, t) is statistically 
homogeneous and isotropic, we obtain 

viJ w - m ' (46) 

where ^{h) is the functional defined in ( |36l) . Thus maximizing If(/i) by solving the Euler- 
Lagrange equation (|371) for /i does double duty: we obtain both a lower bound on the survival 
exponent s and on the ratio r\ (h 2 P) / (h 2 ) . Noting that 

max(/i 2 ) (P) > (h 2 P) , (47) 
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FIG. 9: The dashed curves are the basic lower (|43p and upper (|56p bounds: max{0, (7)} < r/S < 
(v (7 2 ) + (t))/2- The solid curve is the optimized upper bound (|55|) and the dash-dot curve is the 
lower bound (|48p . The circles are the biomass obtained from simulations. The lower bound (|48|) is 
tighter than the basic bound (|4"3"j) only for T < 0.075. 



and using the normalization (h 2 ) = 1, this lower bound also provides a lower bound on the 
biomass !B: 









max 


fa) 



Unfortunately the inequality in (14"T1) is crude, and so the lower bound (T4"8"l) is not always 
tighter than the basic lower bound in (1431) . This is shown in Figure M where (1481) is an 
improvement over (03]) only for T < 0.075. 



B. An upper bound on the biomass 

Having found lower bounds on the biomass in (T4"3"l) and (1481) we now seek a complementary 
upper bound. The first step is to obtain a constraint by multiplying ([1]) by an arbitrary 
positive function f(x) and averaging: 

( V fP 2 -gP) = 0, (49) 

where 

g(x,t) = u- V/ + kV 2 / + / 7 - (50) 
Notice that g{x,t) inherits time dependence from u(x,t) and ^(x, t). 
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We add (3 times ( )49|) onto the definition of 23 

H = (P) + f3( v fP 2 -gP) , (51) 

and then complete the square: 

* (^) ■ 

In passing from the first to the second line we assume that (3 > so that dropping the 
final term in (|52|) results in an upper bound on 23. Minimizing the right-hand side of the 
inequality (153j) yields the tightest bound. The optimal value of f3 is 



P* = \l(l)( g 4)> (54) 



.// \/ 

which is non-negative and therefore consistent with (|53|) being an upper bound. Substituting 
/3* into ( 1531) yields 

For the simple choice f(x) = 1, which implies g(x, t) = j(x, t), (|55p immediately delivers 
the upper bound 

g< W + ^>. (56) 

The upper bound above is sharp in the special case of constant 7 where the stable attracting 
state is 23 = 7/77 and = j 2 /rj. 

To improve on fl56l) requires a better comparison function than /(cc) = 1. One can 
attempt to optimize the choice of f(x) by maximizing the the right hand side of (1551) 
using variational calculus. Unfortunately the resulting Euler-Lagrange equation is very 
complicated and so we compromise by using a simple trial function, such as 

/ = e pcosfciy . (57) 

The adjustable parameter p is determined by minimizing the right-hand side of ([55]) . Since 
p = corresponds to / = 1 this procedure can only improve on (|56|) . This trial function 
procedure has been implemented numerically using the sinusoidal 7(1/) in (J2J) to obtain the 
upper bound indicated by the solid line in Figure [91 In Figure [9] this bound is very tight 
and is the bound closest to the simulation results. 
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FIG. 10: (a) A time series of P(t) compared with two upper bounds; the dashed line is the simple 
upper bound in (|58p and the straight solid line is the optimized upper bound in (|59p using the trial 
function / = exp [0.3960 cos (kiy)]. The dotted line shows the average of the biomass from t* = 100 
to 400. (b)-(c) The open circles indicate the biomass obtained from numerical simulation for the 
indicated parameter values. The dashed lines are the upper and lower bounds in (|58p and the solid 
lines are the bounds in (|59|) with the upper bound p-optimized as in (|57p . The lower bound in (I59p 
only appears in (c) where it is not very tight but does forbid extinction for small 



C. Summary of the bounds 



Our simplest bounds on the biomass are 



1 r 



max{0, ( 7 )}<^<-[V\?) + <7> 



We also derived two bounds involving k*, £7* and the details of 7: 

5 \h 



< B < — 
r] max ( h 2 \ ^ 




(58) 



(59) 
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where h is the solution to ( l37j) and g is defined in ( |50l) . The bounds in ( 1591) can be tighter 
than the bounds in (158|) . but the bounds in ([58]) have the advantage that they only require 
knowing the mean and variance of 7, while the more complicated bounds require knowing 
7 everywhere. 

The productivity is subject to simpler constraints: 



We evaluate and illustrate the bounds ( [58]) and ( 1591) in Figure Panel (a) shows 
an example of a time series of P(t) in an adverse average environment (r = —1); the 
optimized upper bound in ( |59|) . indicated by the solid curve, is obtained by optimizing the 
trial function in fl57|) over p. With = 2, a largish value, even the optimized upper bound 
( 159]) is generously large and the lower bound in (159|) is useless — the best lower bound in 
this case is simply < S. 

Figure fTUlfb) shows that both of the bounds in ([58]) and the optimized upper bound in 
( 159]) become tight as T — > 1 and the growth rate becomes spatially uniform. However, the 
growth rate does not need to be very uniform for the bounds to work well. For example, 
with T = 0.5 the growth rate 7(1/) varies from to 7 max , but despite this large spatial 
inhomogeneity the lower bound in (158]) and the upper bound in f]59]) constrain the biomass 
to lie between 0.57 max /?7 and 0.55087 max /?7. 

The effects of stirring on the biomass are shown in fTOT c): increasing U* loosens the bounds 
and sufficiently large U* drives the plankton to extinction. The optimized upper bound in 
( [59]) is tight provided that U* is not too large and the lower bound in (159|) only appears for 
very small [/*. 

Figure [T07 d) shows that there is a value of (roughly m 10 _1 ) at which the biomass is 
maximized. Unfortunately the optimized upper bound in ( |59l) shows dependence on only 
when diffusion is rather strong and the lower bound in ( |59l) does not appear and so neither 
of the bounds containing are delicate enough to indicate the existence of a maximum 
biomass at a particular value of ac*. 

Finally, the overall the performance of the bounds in Figure [10] is worse than in Figure [9] 
because the simulations in Figure [10] use much larger values of U*. The main conclusion is 
that our bounds on B are tight if £7* < 0(1) and A7/7 < 0(1). In the next section we use 
( 1401) and ( ]42l) to obtain joint inequalities constraining 23 and 7\ 




(60) 
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V. SIMULTANEOUS BOUNDS ON S AND T 



The only information concerning 7 provided by the bounds in the previous section is (I60I) . 
In this section we supplement ( 15 8|) and ([6"Uj) by deriving bounds involving S and T together. 
The first simultaneous bound is obtained by observing that (jP) < 7max (-P), where 7 max is 
the global maximum of ^(x, t). Therefore 

y < 7maxS . (61) 

We obtain a more elaborate simultaneous bound by adding two versions of zero to the 
definition of IP: 

7 = ( 7 P> + a [S - <P>] + /3 [< 7 P> - 77 <P 2 >] , 

7 + /?7 - a 



aT>- fir} 



P 



m ((7 + /j7 _ a)2) 

4/3?7 v ; 



2/377 

Above a and /? are constants. 

If /3 > 0, then we obtain an upper bound on the productivity by dropping the squared term 
containing P from the right-hand side fl62l) . Optimizing this upper bound by minimizing 
over both a and (3 (Appendix [C]) gives 



r£ « s+ C + w-( B -f) 2 ' (63) 

where 

a 2 = < 7 2 > - <7> 2 • 

Returning to fl62|) and taking /3 < we obtain a complementary lower bound by again 
dropping the squared term containing P. Maximizing over a and (3 gives an expression 
analogous to ( 1631) . except that the inequality and the sign of the square root are reversed. 
Thus we have proved that 

P-CB) < 7 < T+OB), (64) 
where the functions y_(S) and y+(S) are 



2?7 V 4?7 2 V 2 V 



y ± (S) s < 7 )S + _±^/iLZ-(S-^) . (65) 



The inequality above constrains the system to fall within the intersection of the first quadrant 
of the (23, y)-plane and the ellipse defined by the arcs CP+(CB) and T_(I>) — see Figure [TIT a). 
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FIG. 11: (a) Bounding ellipses obtained from (j64f) and the sinusoidal model in ([2]); units are 
non-dimensional with 7 max = rj = 1. As T decreases the growth rate becomes more spatially 
inhomogeneous and the ellipse expands and the simultaneous bounds become less tight. However, 
because (P < 7max^> much of the expanded ellipse is inaccessible, (b) The shaded region indicates 
the allowed region if (7) > 0. The point (B, ?) = ^((7) , (t) 2 ) is indicated by A. (c) The shaded 
region indicates the allowed region if (7) < 0. In this case the point A lies outside the first quadrant 
and so the constraint (1681) is toothless. 



Introducing the quantities 



yM±M, s . E VazM, (66) 

2r) 2r) 



and noting from (I56p that 23 < 23 + , we rewrite (1651) as 



2 



y±(S)=r ? [ v /S + (S_ + S)± v /S_(S + -S)j . (67) 

Taking y_(S), the right hand side of (16 7p has a double root at 23 = 0. Thus the 13-axis is 
tangent to the arc IP = T_(23) at the origin of the [25, Jj-plane (e.g, see Figure [TT]) . 
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In ( 1551) and (1601) we obtained the lower bounds 

max {0, (7) } < rf£> and (7) x max {0, (7) } < r/CP . (68) 

I — I 

Employing (165|) . and the upper bound CP < 7 max CB, in concert with (lOTl) 20 we restrict the 
system to the shaded region in Figures ITlT b) and (c). Specifically, the system must fall 
within the intersection of: 

1. the first quadrant of the (23, CP)-plane; 

2. the interior of the bounding ellipse; 

3. the wedge < CP < 7 max CB; 

4. the lower bounds in (165]) . which define a quadrant with southwest apex A in Figure 

nib). 

These joint constraints — which use only (7) , (7 2 ) , 7 max and 77 — provide basic informa- 
tion about the possible range of the biomass and productivity. 

Figure [121 illustrates the joint bound. In panel (a) the parameter U* is varied by a factor 
of 500 while t u = Uk m r in (J5J is fixed. Increasing U* decreases both the biomass and the 
productivity so that for large £/* the system is very near the lower arc of the bounding ellipse 
CP_ in (165]) . Panel (b) shows the effect of varying k*. T < in this case and so large 
drives the system to the origin (extinction). Small also causes the system to head towards 
extinction when [/* = 2. For U* = the system heads to the local solution P(y) = j(y)/T] as 

— > 0. For both values of U* the maximum biomass is achieved with a moderate value of 
/«*: about 0.04 and 0.1 for Z7* = and 2, respectively. The maximum values of the biomass 
corresponding to these parameter values are about 0.23 and 0.19. 



VI. CONCLUSIONS AND DISCUSSION 



Previous work on the FKPP equation and its relatives focused on the survival 



transition 



11 
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231 ] . filamentation 



and front propagation 



QflB 



-extinction 



4j. Here we 



consider the survival-extinction transition and find the existence of a survival region in 
diffusivity-velocity parameter space even if the average growth rate is negative. After con- 
sidering the survival-extinction transition we move on to estimating the average and variance 
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FIG. 12: (a) The joint bounds illustrated with variable U* and fixed r u . £/* increases from 0.1 to 
50 (b) The joint bound illustrated with variable k*; the maximum value of is 1 for both U* = 
and 2. 

of the concentration of plankton (chemical products) in the domain by deriving upper and 
lower bounds on the biomass and productivity. 

Our bounds on the biomass and productivity apply to any plankton model (or chemical 
reaction) with a linear growth term, quadratic damping, and periodic or no-flux boundary 
conditions. The simplest bounds in (1581) do not require detailed knowledge of the growth 
rate, only the mean and variance. Furthermore, the bounds in (1581) do not use information 
about diffusion or the velocity field. This may be attractive, depending on whether or not 
the information is available. 

If the flow is statistically homogeneous and isotropic and (|w| 2 ), the diffusivity and the 
details of the plankton growth rate are known, then we can find more elaborate and possibly 
tighter bounds on the biomass. The upper bound obtained via the trial function procedure 
fl59|) is always tighter than the upper bound in fl58|) while the lower bound in (fl"8l) may be 
tighter than that in (15*81) . Additionally, even when the bound ("""TBI) is not tight, it may forbid 
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extinction, as in Figures 191 and [TOT c) . 

In addition to constraining the biomass, we also derived bounds on the productivity. Our 
definition of productivity, 7 = (jP), is unusual as it includes regions where the growth rate 
7 is negative. Regardless, it is still useful to bound this quantity because from (|40p . (jP) 
equals the mean-squared plankton concentration and may be combined with the biomass to 
find the variance of the plankton concentration. 

Finally, we derived simultaneous bounds on the biomass and productivity which con- 
strain the system to lie inside a certain portion of the biomass-productivity plane. The 
system is close to the boundary of this region for moderate to large diffusivity and velocity. 
These bounds may be useful as an alternative to parameterizing sub-grid scale processes in 
ecological model and for predicting the results of chemical reactions. 

In concluding, we note that the filamentation transition discussed by Neufeld, Lopez 
and Haynes ?J appears in our model if we take T > and make the simplification k = 
0. Examining this transition allows us to assess the effect of small-scale structure on the 
productivity (the biomas is constrained to equal T by ( 1421) ) .The filamentation transition 
occurs when the flow Lyapunov exponent, A in ([7]), exceeds the rate of damping back to 
local equilibrium and P develops narrow filaments which are not smooth in the direction 
orthogonal to the filaments. If the damping back to local equilibrium is strong enough then 
no filaments form and the field remains smooth even if k = 0. Figure [131(a) shows transects 
along the line x = for two different runs using the same sequence of phases <j) x and 4> y in 
([3]) and the same parameter values, except for U* and r*. The argument of the Lyapunov 
function r u in ([8]) is held fixed and so the Lyapunov exponent of the flow increases linearly 
with U*. For the case A = 0.16 the plankton distribution is smooth, but for A = 1.03 the 
distribution is filamented with many sharp peaks in the transect, indicating that we have 
passed throught the filamentation transition. 

Figure fl3T b) shows the Holder exponent a, which is defined in j3] as 

5P~\5x\ a . (69) 

In our model a ~ 1 for A < 0.25 and then smoothly decreases to about 0.4 for A ~ 1.6. In 
between A = 0.25 and A ~ 1.6 the system transitions from smooth to filamented with the 
inflection point of a(A) occurring at A = T, as shown by the dashed line in Figure [131( b) . 
Figure [TBTc) is the same as panel (a) except is small but non-zero. The effect on 
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FIG. 13: (a) Transects along the line x = for two solutions of (JTJ) with k, = using the same 
sequence of phases 4> x and 4> y in ([3]) but different values of A in (J7|). A is varied by holding t u and 
m fixed and varying [/*. The transect with A = 0.16 is smooth, but the transect with A = 1.03 
is not. (b) The Holder exponent (|69p . The vertical dashed line is the line A = T and coincides 
with the inflection point of a (A), (c) Transects along the line x = for two solutions of (pQ) with 
= 10 -4 . The A = 0.16 transect is very similar to the corresponding transect in (a), but the 
A = 1.03 transect is noticeably smoother, (d) The productivity T as a function of A. 

the previously smooth case (A = 0.16) is imperceptible, but the previously filamented case 
(A = 1.03) is noticeably smoother. Panel (d) shows that the productivity decreases as A 
increases . The biomass £ is always 0.75 for = by ( H2l) and constrained by ( 1581) to lie 
between 0.75 and 0.7603 for = 10~ 4 . For the values of A used in Figure [TBI S varies by less 
than 1%. Therefore, the variance, given by r/ _1 — S 2 , decreases as the Lyapunov exponent 
increases and the field becomes less smooth. This counterintuitive result occurs because 
with increased stirring the plankton feel the average growth rate and carrying capacity and 
consequently the magnitudes of the deviations from the mean are decreased as shown in 
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[TBTa) and (c). The filamentation transition has no effect on the biomass or the productivity 
in our model and this may be seen by comparing panels (b) and (d). In panel (b) the 
Holder exponent remains approximately 1 until A ~ 0.25, while the productivity in panel 
(d) begins dropping immediately. This insensitivity to small scale structure indicates why 
crude inequalities, like the ones developed here, may be able to successfully constrain gross 
statistics such as B and 7. 
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APPENDIX A: THE LYAPUNOV EXPONENT 

Under the action of the velocity field ([3]), a particle starting at (x, y) at t — moves to the 
point (V, y') at t = r. This area-preserving transformation can be written down explicitly, 
and one can then calculate the Jacobian matrix which has the form: 

• [ ) ( ■ V (ad 

where (3 = Uk m r/\/2 and s n = sin^, with 9 n a random phase. We denote the transpose of 
J by JT. 

Consider an ensemble of infinitesimal line elements; at t — each element has length £ , 
and is oriented along a randomly directed unit vector e. After one iteration, an element is 



stretched by a random factor e T J 1 r J 1 e to a length £ 1 = e 1 Jj Jie£ . After n iterations the 
length is given by 



4 = ^/eT n ---JlJ 1 ---J n e£ . (A2) 

The matrix Sn = ( J\ ' ■ ' Jn) 1 (Ji ■ ■ • Jn) is real and symmetric with determinant one. Hence 
we have the representation 

d n = Q J EQ , (A3) 
where the matrix Q corresponds to rotation through a random angle \ and 

(A4) 
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a > 1 is an eigenvalue of Sn- The Lyapunov exponent A defined in ([7]) is given by A = 
(ln(£ n /£o)) / (nr), where () indicates an average over the random angles in e and Q and over 
the the random eigenvalue a. We cannot calculate this average analytically because the 
average over a is complicated. However the average over the angles is a standard integral: 

(ln(4/4)) x = \ (In {a cos 2 X + oc~ l sin 2 X ) ) x , (A5) 
= In [± (a 1 ' 2 + CT 1 / 2 )] (A6) 

An efficient Monte Carlo procedure takes advantage of (1A6I) by generating many realizations 
of 3 n , calculating a for each realization and obtaining the Lyapunov exponent as 

A = ^<ln[H« 1/2 + «- 1/2 )]>- (A7) 

This procedure using 4 x 10 4 realizations of Sn with n = 20 gives the solid curve in Figure 
Ha). 



APPENDIX B: APPROXIMATE SOLUTIONS OF A MATHIEU EIGENPROB- 
LEM 

In section 3 we solve a variety of eigenproblems related to Mathieu's equation, written 
here on the domain —it < y < it as 

KQ yy + (cosy - E)Q = . (Bl) 

This is not the standard form of Mathieu's equation, but it is convenient for the problems we 

face in section 3. Requiring that the solution of (1B1I) be periodic determines an eigenrelation 

K(E) between the parameters K and E. In this appendix we focus on the first mode (that 

is the ground state) and obtain the main features of the function K(E) via perturbation 

theory pivoted round the complementary limits K — > and K — > oo. This result is used to 

deduce the simple analytic approximations indicated by the dashed curves in Figure [5fb). 

If K ^> 1 we define e = K^ 1 1 and regard E{e) as an eigenvalue. We then solve (IBlj) 

via a regular perturbation expansion in e. The result is 

e 2 e 
Q = l + ecos?/ + -cos2y + 0( e 3 ), E = - + 0(e 3 ) , (B2) 
o 2 

or 

K ~ ^ , ifK»l. (B3) 
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If K <C 1 then E is close to 1 and it is convenient to define a parameter 5 <C 1 by 
E — 1 — 5 2 . Thus in flBlj) the growth rate cosy — E is positive only in a small region of size 
S centered on y = 0, and within this oasis the equation is 

KQ yy + (5 2 - y 2 /2) Q = 0(6 4 ) . (B4) 

The solution of this quantum oscillator problem is standard: Q ~ exp(— y 2 /A5 2 ) and K pa 
25 2 , or 

K pa 2(1 — E) 2 , iftf<l. (B5) 
The two approximations in (IB3I) and (IB5I) can be combined into a single expression 

Km Hx±*S£*L. m 

The above is asymptotically exact as E —>■ and as E — > 1, and interpolates the function 
K(E) over the range < E < 1 with maximum error of about 4.5% at E pa .67. 

The first application of (IB6I) is to the eigenproblem obtained by setting = in (T23l) . 
Making the identifications if = «*/(l — T) and = —T/(T — 1) we obtain ( 1B1I) . Then 
rewriting (1B6|) in terms of and |T| = — r we obtain the simple approximation (|24[) . 
Figure \50Jd) compares this approximation with the numerical solution obtained by adapting 



program 21 of Trefethen 25]. 



APPENDIX C: OPTIMIZATION 



We continue from f[62|) filling in some algebraic details. Dropping the square containing 
P from the right-hand side yields: 

y<^+ ((7 + ^- a)2) - (ci) 

Minimizing the right hand side of (IClj) over (3 we find the optimal value of (3: 



Plugging /3* into fICll) yields 



3 , <nfl+ (7(7-a)>-t-V(7»X(7-a)') - (C3) 

2r7 
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(Note that /?* is positive and so our earlier assumption that [3 > is valid and we have found 
an upper bound.) Next we minimize the right-hand side of flC3j) over a. After some work, 
we find that the optimal a* satisfies 

(7) - 



where 



A , (C4) 



A.^Us-M, (C5) 



The inequality (1351) ensures that < .A < 1. 

Solving flC4p . and taking the branch consistent with A > 0, we find that 

Aa 

a * = w -7CT (C6) 

where 



^V<7 2 >-<7> 2 . (C7) 
Substituting (1C6P into (103|) we obtain 



2 / / 9\ / / \ \ 2 



W _ U - M 

2?7 ' " y 4?7 2 V 2r ? 



y<(7)S + - + a t /^-(S-^) . (C8) 



Repeating the above procedure with (3 < gives the other half of the bounding ellipse. 
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